arXiv:1505.03232vl [physics.optics] 13 May 2015 


SEAMLESS INTEGRATION OF GLOBAL DIRICHLET-TO-NEUMANN 
BOUNDARY CONDITION AND SPECTRAL ELEMENTS FOR 
TRANSFORMATION ELECTROMAGNETICS 

ZHIGUO YANG * 1 , LI-LIAN WANG 1 , ZHIJIAN RONG 2 , BO WANG 3 AND BAILE ZHANG 4 


Abstract. In this paper, we present an efficient spectral-element method (SEM) for solving 
general two-dimensional Helmholtz equations in anisotropic media, with particular applications in 
accurate simulation of polygonal invisibility cloaks, concentrators and circular rotators arisen from 
the field of transformation electromagnetics (TE). In practice, we adopt a transparent boundary 
condition (TBC) characterized by the Dirichlet-to-Neumann (DtN) map to reduce wave propaga¬ 
tion in an unbounded domain to a bounded domain. We then introduce a semi-analytic technique 
to integrate the global TBC with local curvilinear elements seamlessly, which is accomplished by 
using a novel elemental mapping and analytic formulas for evaluating global Fourier coefficients 
on spectral-element grids exactly. 

From the perspective of TE, an invisibility cloak is devised by a singular coordinate transfor¬ 
mation of Maxwell’s equations that leads to anisotropic materials coating the cloaked region to 
render any object inside invisible to observers outside. An important issue resides in the imposi¬ 
tion of appropriate conditions at the outer boundary of the cloaked region, i.e., cloaking boundary 
conditions (CBCs), in order to achieve perfect invisibility. Following the spirit of m , we propose 
new CBCs for polygonal invisibility cloaks from the essential “pole” conditions related to singular 
transformations. This allows for the decoupling of the governing equations of inside and outside 
the cloaked regions. With this efficient spectral-element solver at our disposal, we can study the 
interesting phenomena when some defects and lossy or dispersive media are placed in the cloaking 
layer of an ideal polygonal cloak. 


1. Introduction and problem statement 

Accurate simulation of wave propagations in inhomogeneous and anisotropic media plays an 
exceedingly important part in a wide range of applications related to the exploration and design of 
novel materials that enjoy unusual and remarkable properties in steering waves. In many situations 
involving time-harmonic wave propagations, the development of high-order methods (i.e., spectral 
and spectral-element solvers) for the Helmholtz equation and time-harmonic Maxwell equations, is 
of fundamental importance. 
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We are concerned with the two-dimensional Helmholtz equation governing time-harmonic wave 
propagation in anisotropic media: 

V • ( C(r ) V«(r)) + k 2 n(r)u(r) = f(r ), r = x = (x,y) € R 2 , (1.1) 

where k > 0 is the wave number in free space. In general, we make the following assumptions. 

(i) C is a symmetric positive definite matrix in R 2x2 , and for some positive constants co,ci, 

0<Cb<£‘C£<ci, V£ G R 2 , a.e. in l 2 . (1.2) 

(ii) The coefficient 

0 <n< n \, a.e. in R 2 . (1.3) 

(iii) The inhomogeneity of the medium is confined in a bounded domain with Lipschitz 
boundary, and / is compactly supported in disc B R of radius R > 0 (see Figure [P] ): 

C = I 2 , n= 1 in R 2 \S2_; supp(/) C B R , (1.4) 

where J 2 is the 2x2 identity matrix. In what follows, we are interested in the case where 

We impose the well-known Sommerfeld radiation 
boundary condition upon the scattering wave: u sc := 
u — U[ n (where u- ln is a given incident wave): 

d r u sc — \k u sc = o{r~ 1 ^ 2 ) as r —> oo, (1.5) 
where i = y/—l is the complex unit. 

The challenges of the above problem are at least 
threefold: (i) unboundedness of the computational do¬ 
main; (ii) indefiniteness of the variational formulation; 
and (iii) highly oscillatory solution decaying slowly 
when t > 1. In addition, the coefficients C(r) and 
n(r ) might be singular at some interior interface in 
S2_ (see Section |3j). 

The methods of choice to deal with the first issue typically include the perfectly matched layer 
(PML) technique [5], boundary integral method [221 [30] , and the artificial boundary condition 
HH1E3 E3 HU- The latter is known as the absorbing boundary condition (ABC), if it leads to 
a well-posed initial-boundary value problem (IBVP) and some “energy” can be absorbed at the 
boundary. In particular, if the solution of the reduced problem coincides with that of the original 
problem, then the related ABC is dubbed as a transparent (or nonreflecting) boundary condition 
(TBC) (or NRBC). In this paper, we adopt the exact TBC (see [33] and Figure 


d r U sc '-^R [^sc] — 0 &t 

rB, 

(1.6) 

where the DtN map is defined as 



OO 

^rW\ = X! with 


(1.7) 

- i r 27T 

^m=^J o 1p(R,9)e~ ime de, Tm 

kH$\kR) 

H { rk\kR) ■ 

(1.8) 


Here, Hm\z) is the Hankel function of the first kind (cf. [T]). This yields the exact boundary 
condition for the total field: 


d r u - &r[u] = d r u iB - ^R[uin\ := h at T fi . 

We find it is advantageous to impose DtN TBC for the following reasons. 


is a penetrable scatterer. 


DtN 



i 

i 

i 

'T, 


/ x R 


FIGURE 1.1. Illustration of geometry 


(1.9) 
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(i) The original problem in R 2 reduces to an equivalent boundary value problem (BVP) in Br. 
One can place as close as possible to S2_, as long as the inhomogeneity of the media and 
support of the source term are confined in Br. 

(ii) It is essential for accurate and stable simulations especially when the wavenumber is large. 


However, the TBC (1.9) is global in space, that is, evaluating &r[u](x) at any point xq € Tr 
requires to compute the global Fourier integral along the circle T#. This poses challenges in solving 
the reduced Helmholtz problem by a local element-based method. 

(i) The spectral-element solution (using curvilinear elements along T/j) is piecewise continuous 
(i.e., only in C^T^)), and defined on spectral-element grids, so the interplay between Fourier 
points and spectral-element grids via interpolation and fast Fourier transform (FFT) only 
leads to a first-order convergence. 

One can evaluate the Fourier integral by a composite rule with a decomposition coherent 
to the spectral-element partition, but due to the elemental mapping between the curvi¬ 
linear element and reference square, a numerical quadrature is usually necessary, which is 
prohibitive as the integrands are highly oscillatory for high Fourier modes (see ( 2.28| )). 

One main purpose of this paper is to seamlessly integrate the global DtN BC with local spectral 
elements. The key idea is to construct a new elemental mapping between the curvilinear elements 


(ii) 


along r# and the reference square (see Figure 2.1), which leads to exact evaluation of the Fourier 
integrals. It is noteworthy that Fournier [2] proposed a method for calculating global Fourier 
coefficients for given nodal values on non-conforming spectral elements, where a similar semi-analytic 
approach was essential for the success of the method therein. We also remark that the recent work 
[TO] addressed the integration of one-dimensional DtN TBC imposed on a line segment with standard 
rectangular elements along the boundary. Different from these works, our “local-to-global” method 
is built upon the use of curvilinear elements seamlessly fitting the circular boundary, and the design 
of a new elemental mapping leading to exact calculation of the involved Fourier integrals (see 
Subsection 12.3.21). 


Underpinned by the advent of metamaterials, transformation electromagnetics (TE) (cf. [361 
m) provides a powerful tool for creating novel devices and new materials with unconventional 
properties (see, e.g., [Ml [Ml HI HOI 113 IM| and [44j for many original references therein). Some 
exciting applications of TE include the invisibility cloaks (see, e.g., [Ml EH) , rotators (see, e.g., 0) 
and concentrators (see, e.g., [39]), which naturally give rise to the model problem (|l.l|)-( 1.5). In 


particular, the invisibility cloak is devised by a singular coordinate transformation [36] that leads 
to singular materials coating the cloaked regions and preventing waves from penetrating into the 
inside region. The imposition of appropriate interface conditions at the inner boundary, i.e., CBCs, 
where the material parameters are singular, becomes critical. Significant efforts have been devoted 
to CBCs for circular cylindrical and spherical cloaks. Ruan et al. m first analytically studied the 
sensitivity of the ideal circular cloak |36| to a small (5-perturbation of the inner boundary. Zhang et 
al. [M] provided deep insights into the physical effects of the singular transformation (also see 52) ■ 
To shield the incoming waves, the perfect magnetic conductor (PMC) condition was imposed at the 
inner boundary in finite-element simulations (see, e.g., [Til 1251151] ). Weder [13] proposed CBCs for 
the ideal spherical cloak of Pendry et al. [M] from the perspective of energy conservation. Lassas 
and Zhou [26j[25] proposed some non-local pesudo-differential CBCs. Based upon the principle that 
a well-behaved electromagnetic field in the original space must be well-behaved in the transformed 
space as well, Yang and Wang [31] obtained CBCs for circular and elliptical cloaks that intrinsically 
relate to the essential “pole” conditions of a singular transformation. 

The polygonal cloaks enjoy more flexibility to hide objects with complex shapes, which are 
however much less studied. Indeed, many of the previous principles and approaches for CBCs 
are not extendable to the polygonal case. Following the spirit of [48], we propose new CBCs under a 
“local” coordinate system (see Proposition 3.1), under which the governing equation in the cloaked 


region is decoupled from the exterior region. Accordingly, no wave can propagate into the cloaked 
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region, and vice versa. We emphasise that the new CBCs are indispensable for spectrally accurate 
simulations. We also show that the proposed spectral-element solver provides a reliable tool to study 


how the defects affect the perfectness of an ideal cloak (see Subsection 3.5). 


The rest of the paper is organised as follows. In Section[2j we review the form invariant of Maxwell 
equations and illustrate the derivation of the above model problem. We then introduce the technique 
to seamlessly integrate the global DtN BC with local spectral elements. Section [3] is for accurate 
simulation of polygonal invisibility cloak, where new CBCs are derived and efficient techniques 
are introduced to deal with singular material parameters. Various numerical results are provided 
to show the perfectness of invisibility, and the effects of defects and lossy or dispersive media. 
Section [4] concerns the extension of the spectral-element solver to the simulation of electromagnetic 
concentrators and rotators. 


2. TE AND SPECTRAL-ELEMENT DISCRETIZATION OF DtN BC 

In this section, we first illustrate the scenarios of the aforementioned model Helmholtz problem 
arisen from transformation electromagnetics. We then discretise the model problem by spectral- 
element method and focus on how to seamlessly integrate the global DtN BC with local elements. 

2.1. Form invariant of Maxwell equations. Consider the time-harmonic Maxwell system: 

Vf x E — icD/xo H = 0, Vf. x H + iweo E = 0, (2-1) 

in Cartesian coordinates: r = x = ( x,y,z ) £ M 3 , where the electric permittivity eo, the mag¬ 
netic permeability y o, and the angular frequency a; are positive constants. Note that e~ lut time- 
dependence is assumed for the electric and magnetic fields. 

A remarkable property of the Maxwell system is its form invariant under any coordinate trans¬ 
formation (cf. [3Z]). More precisely, given a coordinate transformation r = r(r) with the Jacobian 
matrix J = dr /dr , the transformed Maxwell system takes the same form: 

V x E — iw/Lto/x H = 0, V x H + iweoe E = 0, (2-2) 

where V x is the curl operator in the new coordinates, and 

E(r) = (J t )~ 1 E(r), H (r) = (J*) -1 H (r), fi = e = JJ ( /det(J). (2.3) 

We are concerned with the two-dimensional electromagnetic wave propagations in media with 
in-plane anisotropy. Accordingly, under the transverse-electric (TE) polarization, we consider E = 
(0,0 ,u(x,y)Y and H = (Hi, H- 2 , 0) ( . Letting z = z in the coordinate transformation, the material 
parameters in (2.3) reduce to 


fi = e = 


c 

0 * 


0 

n 



On 

C 12 

0 


C\2 

C 22 

0 


where 


C = 


T T 

J ai 4 cn 


1 


det(Jcn) ’ 
Note that det(C) = 1, and 


n = 


det(J cn ) 


with Jen := 


d$x djjX 
dxV dfjy 


fi 1 = e 1 = 


Then we derive from the first equation of 


H = 




-1 


-1 


C*22 

-Cl 2 

0 

and 
1 


-Ci 2 

Cn 

0 


(2.4) 


(2.5) 


( 2 . 6 ) 


that 


V X E . (^ 2/5 . (C'i2'U / x H - C'll'U'x ^'12^'yT^) • (^•^0 

lCJ/iQ !Cc7yL/-o lCJ/iQ 
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Inserting it into the second equation of (2.2 1 , we obtain the two-dimensional Helmholtz equation: 

V • (C(r) Vu(r)) + k 2 n(r ) u(r) = 0, (2.8) 

where k = ujy/e o/io is the wavenumber in free space. 

In Sect ions [3]4j we shall introduce the coordinate transformations for polygonal invisibility cloaks, 
concentrators and rotators, and compute the corresponding material parameters C and n via ( |2.5| ). 
It is noteworthy that in all cases, the coordinate transformations are identity in R 2 \ (f2_ U S2 + ) 
(cf. Figure B so ( |1.4[ ) can be met. Moreover, we can derive ( |2.12| below from the standard 
transmission conditions (see, e.g., [35] Sec. 1.5] and (32]), that is, the continuity of the tangential 
components of E and H at the interface F := <9H_. 

In summary, the problem of interest reads 


V • (C(r)Vu(r)) + k 2 n(r)u(r) = f(r) in Br, 
[w] = [CVti] = 0 at r, 
d r u — £?r[u] = h at T/j, 


where 


m := u — u 


[CVu] := n ■ (C~Vu~ - C + Vu+), 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 


:= u\n ± , C ± := C|n ± and n is the unit outer normal vector along T. 


2.2. Spectral-element scheme. Let H be a generic bounded domain, and L 2 (H) be the space of 
square integrable functions with the inner product and norm denoted by (•, -)n and || • ||n as usual. 
The Sobolev space H m (£i) with to > 0 is defined as in Admas [2J with the normal ||u||m,n- Define 
the trace integral 

(2.13) 


{u,v)r R ~f uvd'y. 
J r R 


A weak formulation of ( |2.9[ )-( 2.11) is to find u € H 1 (Br) such that 

&(u,v) : = (CVu, S7v) Br - k 2 (nu,v) BR - {&r[u],v)t r 

= &(v) := -(f,v)s R + (h,v)r R , Wv £ H\B r ), 


(2.14) 


where by (1.7)-( 1 . 81 , 


j-) OO p2,7T p2,7T 

{Mu\,v) Tr = — Y, T m (j o u(R,8)e- ime dd)(J^ v(R,8)e~™°dff). (2.15) 


| m |—0 


Remark 2 . 1 . Recall that (cf. 

ItoI +1 


P. 87]) 


R 


<Re(7^)<—, 0 < Im(T m ) < k, |m|= 0 ,l," 

K 


(2.16) 

□ 


We can claim the unique solvability of (2.14) from ( 1 . 2 )-( [T~4| ) and ( 2.16| ) (cf. [[55] 1. 

For simplicity, we assume that the scatterer is a polygonal domain or a disk, though our 
approach is extendable to more complicated domain. We partition the computational domain Br 
into a finite number of non-overlapping straight-sided or curvilinear quadrilateral elements {kl e }f =11 
such that the inner interfaces are aligned with the “edges” of the elements. In particular, we have 


e r 


= U T n = U 0 e+i] = [°. M r 

e—1 e=l 

where I]f := n ^ 0 for all e € {l,-- - ,E}, and 9i = 9e r +i (see Figure 
X e : Q ■= (—1, l) 2 — > kl e be a one-to-one elemental mapping defined by 

r = x = (x,y)= X e (C,r?) := (x?(£>»?),* 2 (£,»?)). V(^,??) G Q. 


2.1 


(2.17) 
(a)). Let 


(2.18) 
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Recall that one commonly-used elemental mapping, originally proposed by Gordon and Hall HE 
transforms Q to any quadrilateral fl e with straight or curved sides (see, e.g., mm )• Here, we shall 
use a special Gordon and Hall transform in (2.23)-(2.24) below. 


Denote by Vn the set of all polynomials of degree at most N in [—1,1], Introduce the spectral- 
element solution space 


■= {v £ C(B r ) : v(ic)|n* = v(x e ) £ V%, 1 <e<E). 
The spectral-element approximation of (2.14) is to find u^ £ V^ such that 


@{u%,v%) = Vv%£Vft. 


(2.19) 


( 2 . 20 ) 


In view of Remark 2.1 we can show the well-posedness of (2.20) as with (2.14). 


2.3. Seamless integration of SEM with DtN TBC. As usual, the continuous inner product 
(•,-)n can be evaluated by element-wise discrete inner product based on tensorial Legendre-Gauss- 
Lobatto (LGL) quadrature, and likewise for the term (h, Vn)t r (see, e.g., [T2]). However, much care 
is needed to deal with the term {&r[un\, Vn)t r , as the DtN operator is global, but the spectral- 
element solution is piecewise. One can evaluate (2.27) by using the fast Fourier transform (FFT), but 
this requires an intermediate interpolation to interplay between spectral-element grids and Fourier 
points. Since uf/\r R £ C° , a naive interpolation only results in a first-order convergence. 

In what follows, we introduce an efficient semi-analytical means to compute {^r\u%\,Vn)t r - Let 
{£j = r/j}jL 0 (i n ascending order) be the LGL points in [—1,1], and let {lj}jL 0 be the associated 
Lagrange interpolating basis polynomials. Correspondingly, the spectral-element grids and basis on 
Cl e are given by 

Xij =X e (&,Vj), =h{Olj(v), 0 <i,j<N. (2.21) 

Formally, we can write 


Un(x,v) | ne =J2u e ijh(Olj(v), 


( 2 . 22 ) 


where the unknowns {u® } are determined by the scheme (2.20). 



FIGURE 2.1. Curvilinear elements and tensorial LGL points on the reference square and 
a curvilinear element via the new elemental mapping. 

We now particularly look at the Gordon-Hall transform for a curvilinear element Q e with vertices 
{(xf,yf)}j = i along F/i (with three straight sides). Let {7T®(t),f £ [— l,l]}j =1 be, respectively, the 
parametric form of four sides such that 

=nl{l), tt?(1) = 7ri(l), tt!(-1) = 7r|(l), tt|(-1) =tt|(-1), 


(2.23) 
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see Figure 2.1 (b). In this case, the Gordon-Hall transform takes the form 
x = X e (C,V) = 7 r l(0^-7T^ + 7r 3 (£)Er^ + 


<(»?) 


- 

- (E(-l) 


1-4 

2 

1-4 




1 + A 1 + V 


(2.24) 


*I(1) : 


l + 4\ 1 — »7 


2 dW 2 

where the edge r] = 1 of Q is mapped to the arc 1 ^ = {r = R, 9 £ (0 e , 9 e+ 1 )} of S7 e , i.e., 
r R : ® = Xi(£,l) =7rn(0» V = X 2^,!) = 7Ti 2 (0) V£e(-l,l). 


Accordingly, the spectral-element grids in polar coordinates on 1]^ (see Figure 2.1) satisfy 
cosdj = i? _ 1 7 r® 1 (^j) or sin#| = i? _1 7r® 2 (Ct)> 1 < j < N. 

We now turn to (^[u®],trj|)r H in (2.20). Thanks to (2.151 and (2.22), we need to evaluate 


(2.25) 

(2.26) 


r 2ir 


Er r Q e 


l N 


0,2/) I 


~ ime de = J2 


l N 


0,2/) I 


„—imQ 


d 6 


e=l ' 
Er 


= EE 


u iN 


U(Z)e 


1-1 


-im0(O^ d £ 

d4 


(2.27) 


As the nodal basis {/;} can be represented in terms of Legendre polynomials, it suffices to compute 
Km ■■= f Pn{0 e~ imm ^ d£, for n > 0 , \m\ > 0 , (2.28) 

where P n is the Legendre polynomial of degree n, and by ( |2.25| ), 

^ ~ = W^frO )] 2 + [O^O )] 2 ■ (2-29) 

It is seen that the integrand is highly oscillatory for large |m|, and the efficiency and accuracy in 
computing Km essentially relies on the choice of the parametric form for 7 Ti( 4 ). We next introduce 
a parametric form that allows for exact evaluation of (2.29) by analytic formulas (see Propositions 
2.1|2.2 ). To stimulate the idea, we first consider a commonly-used parametric form. 


2.3.1. A commonly-used parametric form for (£). Following the ideas of the cubed-sphere trans¬ 
formation (cf. 0011521) and the “ray” coordinates (cf. [23]), one can project the secant line: {x\,yf), 
(x 2 , y 2 ) to the arc T)) via the “rays” from the origin. This leads to the parameterisation: 

Rdi{0 Rd 2 (0 


<(4) = «1 (4)^(4)) = 


where 


di(0 = 


'H + x -P- 


Vdm+dm’ vdm+dKoJ' 


j yl-vl? , 2/1 + 2/! 
“ 2 (4) = —o—4 H —• 


Since cos 0 = I? 1 7 r® 1 ( 4 ), we find 

m = 


fa, 

if 9 € [0,7r), 

-i ( d^) \ 

a := cos —7 ^ 

[271 — a, 

if 9 € [ir, 27r), 

w<% (o+d 2 mJ 


and (2.29) reads 


d9 \xM-x* 2 yt\ 


d4 2(df(0 + ^(0)’ 


V 4 e [-1,1]. 


(2.30) 

(2-31) 

(2.32) 

(2.33) 


Inserting (2.32) and (2.33) into (2.28), one immediately finds that Km appears complicated and 
must be evaluated numerically. However, the integrand is highly oscillatory, when |m| is large. 
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2.3.2. A new parametric form for 7T®(£). We next take a very different route to parameterise rj). 
The essential idea is to look for 


* 1(0 = (*ii( 0 >* 12 ( 0 ) = #(cos 0 ,sin 0 ), 0 €[ 0 e , 0 e+ i], ^ e [— 1 , 1 ], 

such that d j = a d£, that is, the arc length 7 is linear in £. 


(2.34) 


Proposition 2.1. Let Q e be the curvilinear element as in Figure 2.1 (b). Then the new elemental 
mapping from the reference square Q to Ll e takes the form 


where 

with 


*=^(e>^ + 

, = vy.0 1 -^ + + tti* 

r l(0 = (*u(0>*12(0) = R{ COS(0 e £ + Pe), Sm(0 e £ + /3 e )) , 


f. 0e+l — 0e a 0 e + 9 e +l 
vp — ~ 5 Pe — 


(2.35) 

(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 

Inserting it into (2.34) leads to the new parametric form (2.37). Then we obtain (2.35)-(2.36) from 
the equations of the straight sides, e.g., 


2 2 
Proof. Let 7 = at; + b. The arc length along TJ) is 7 = R(6 — 0 e ), so we have 

a£ + b = R(6 — 0 e ). 

Since 0 = 0 e (resp. 0 = 0 e +i) is mapped to £ = — 1 (resp. £ = 1), we find 

a = b = § e R, 0 = 0 e £ + &, £ e [-1,1]. 


<(v) = — 


-V- 


■; V £ [ 1: 1]; 


and the Gordon-Hall transform (2.24). 


(2.41) 

□ 


Observe that in distinctive contrast to (|2.32|)-(|2.33|), the new transformation has a linear depen- 

(2.42) 


dence of 0 in £, so in (2.28), 


d0 


0(0 = #e£ + Pe, = 0 e - 


This leads to the following analytic means for computing the integrals of interest. 


Proposition 2.2. Under the new transformation in Proposition 2.1 the integral in (2.28) can be 
computed by 

Ko = 2 OeSno; Km = ^./^J n+ 1 / 2 (m 0 e ) e~ im ^, (2.43) 

1 y 2mO e 

and K _ m = (Km)* f or n> 0 and m > 1, where J n+ 1/2 is the Bessel function of the first kind, and 
0,/3 e are the same as in (2.38). 


Proof. By (2.281 and (2.42), 


, d 0 


Km = J Pn( £) - d£ = 0 e e~ ir ‘ 


Pn(0e 


—i mQ e £ 


d£. 


(2.44) 


i-i 


It is clear that for m = 0, we have Iqo = 20 e , and by the orthogonality of Legendre polynomials, 
we have I ® 0 = 0 when n > 1. Moreover, we have I® _ m = (Km)*> so we on ly nee d to compute the 
integrals with m > 1. Recall the identity (cf. (3j) 
r 1 


Pn (£) e 


-imx{ jc _ j: 


1 / 2 tt 


i-i 


r L V mx 


Jn+ l/2( mx )i for mx > 0- 


(2.45) 
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Thus, (2.431 follows immediately. 


□ 


2.4. An illustrative numerical example. As a by-product, the new parameterisation provides 
an efficient means to compute the Fourier coefficients via piecewise Legendre approximation. In 
a nutshell, we partition [0, 27r] into {[9 e ,0 e +and approximate the underlying function on 
each subinterval by Legendre polynomials using ( |2.37| so that the analytical formula (2.43) can be 
applied. Here, we provide an example to illustrate this numerical-analytic approach. 

Consider the Fourier expansion of a plane wave (cf. jT) P. 360]): 


ik(x cos OQ+y sinOo) _ 


AmO 


with g m = i m J m (kR)e 


— im#o 


(2.46) 


for some constant 9q, where J, 
-Er 


= n 

| m |—0 

is the Bessel function of the first kind of order m as before. Let 

Note that 

m| increases. We depict the errors against N (with fixed Er = 4) in 

right), 


g m R N be the numerical approximation to g mi and denote the error max| TO i< M |g m — g m R N 
g m decays exponentially as 

Figure [272] (left), and against the number of elements Er (with fixed N = 10) in Figure 2.2 


for R = 1, M = 20, k = 10, 20, 30 and 9 0 = tt/ 4. We observe exponential convergence in both cases. 




FIGURE 2.2. Numerical error maxi m |< A r \g m — g^ R N \ with fixed M = 20, 9o = 7t/4 and 
k = 10, 20, 30. Left: errors against N with Er = 4. Right: errors against Er with N = 10. 


3. Accurate simulation of polygonal invisibility cloaks 


As already mentioned, the invisibility cloak is one of the most exciting examples of transformation 
electromagnetics outlined in Subsection |2.1| In this section, we apply the spectral-element solver 
proposed in Section[2]to simulate the polygonal invisibility cloak, and numerically study the effects of 
defects, lossy media or dispersive media in the cloaking layer. We particularly address the following 
two important issues. 

(i) How to impose appropriate cloaking boundary conditions at the boundary of the cloaked 
region to perfectly hide the objects inside the cloaked region? 

(ii) How to efficiently treat the singular material parameters in spectral-element discretisation 
to ensure accurate simulation? 


3.1. Coordinate transformation and material parameters. As with Pendry et al. [3B], the 
underlying coordinate transformation for a polygonal cloak blows up the origin O in the original 

(a) to the polygonal domain IF = A p B p ■ 


(x, i/)-coordinates in Figure 


3.1 


3.1 


(b), 


• F p in Figure 

which forms the “cloaked region”. It is expected that the waves from outside can not propagate into 
fP so that any object inside is concealed. Accordingly, the polygonal domain f2_ (i.e., the polygon 
AB ■ ■ ■ F) is compressed into the polygonal annulus f2“ = fl_ \ fP , called the “cloaking layer” 
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FIGURE 3.1. Schematic geometry of a polygonal cloak, (a) The polygonal domain in 
the original coordinates (x,y). (b) Through (3.21, the origin is spanned into the polygonal 
domain = A P B P ■ ■ ■ F p that forms the cloaked region. Consequently, the original 
polygonal domain in (a) is compressed into the polygonal annulus f2“ (i.e., the shaded 
part) that forms the cloaking layer, (c) The “local” coordinate system (r,n). 


In fact, the coordinate transformations can be best characterised in two polar coordinates as 
depicted in the diagram: 


^Original: (a:,j/)j 4- > Polar: (r, 0) 


Transform 


-> (Polar: (r,9)j< ->■ ( Physical: (x,y) j (3-1) 


With this, the coordinate transformation for a polygonal cloak takes the form (see, e.g., [53]): 

(r = (1 — p)f + Ri, fe[0, R 2 \, re[R 1 ,R 2 \, 

[0 = 9, 9,0 G [0,2tt), 

where 


(3.2) 


OAp OBp p 

P= -OA = OB="-=OF’ ° <P<1 ' (33) 

and (Ri(9),9) (resp. (R2(6),9)) is the polar parametric form of the side of the cloaked region fF 
(resp. the domain fi_). 

We find it is more convenient to introduce a “local” coordinate system to represent Ri (* = 1,2) 
when we impose the cloaking boundary conditions later on. Consider any side of the insidemost 
polygonal domain fF , say A p B p in Figure 3.1 (c), with vertices and ( x 2 ,y 2 )■ Then its unit 


tangential vector r and unit normal vector n are given by 

(x 2 - xi,y 2 - yi) 


V ( x 2 - Xi) 2 + (y 2 - yi) 2 


■=(n,r 2 ), n = (t 2 , —ti), 


(3.4) 


respectively, which form a “local” coordinate system. Noting that for any (, x , y) = (Ri cos 9, Ri sin 9) 
along this side, we have (x, y) • n = ( x 2 , y 2 ) ■ n = {x\,y\) ■ n, and 


Ri(0) = 


t 2 x 2 - Tiy 2 
t 2 cos 9 — n sin 6 ’ 


R 2 {9)= p-^R^O). 


(3.5) 


We can use (2.51 and (3.2) to derive the coefficients C and n in the Helmholtz equation (2.9). 


More precisely, the coefficients in Hr = f l p _ U f2“ U H + take different forms as follows. 
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(i) In the cloaking layer fi“ , the entries of C in (2.4)-(2.5) are given by (see Appendix |A|) : 

1 / di?i x 


C n — 


Co 2 — 


r — R\ x 2 


r r 2 + r(r — i?i) \ dd r 

r~Ri y 2 + 1 f dRiy 

r r 2 r(r — i?i) \ d dr 


(3.6) 

(3.7) 


and 


C12 — 


r — R[ xy 


1 / di7i x \/'dRiy \ 

r(r — R\) \ dd r / \ d0 r / ’ 


r — R\ 
r(l — p) 2 


(3.8) 

(3.9) 


(ii) In 17+ = Br \ f7_, we have C = I 2 and n = 1. 

(iii) Following [36], we set C = 1 2 and n = 1 in flX 



3.2. Transmission conditions and new CBCs. For notational convenience, we partition the disk 
Br into a finite number of non-overlapping subdomains based on the nature of material parameters. 
Take the setting in Figure [Xl] as an example, and decompose 

B R = fl p _un a _U 17+, 0“ = ufLifi l T , 17+ = Uf =1 17\ (3.10) 

where {17^} are trapezoids. Denote by T+ the “outside” and “inside” boundaries of the cloaked 
region 170, respectively. We further denote by T a := 17+ n 170, i.e., the outer boundary of 17“. 

We next impose interface conditions: (i) at T“ and sides in the radial direction of the trape¬ 
zoids {17^}; and (ii) at the cloaking boundary 1 + . In the former case, the material parameters are 
bounded below and above, so we impose the usual transmission conditions. Recall that under the 
TE polarisation, E = (0, 0, u)‘, so we have 


[u] = [CVit] = 0 for case (i). 


(3.11) 


In the latter case, the material parameters degenerate at T+ (cf. (3.6)-(3.8)), so the above 


transmission conditions are not applicable. In fact, the tangential component of E is not continuous 
across the cloaking boundary. Indeed, it is the singular non-isotropic medium in the cloaking 
layer 17“ that offers the possibility of achieving invisibility within the innermost polygonal domain 
170. In practice, the perfect conduct shell, i.e., PEC (under TM polarisation) or PMC (under TE 
polarisation) was enforced at P+ in simulations (see, e.g., |TT] for circular cylindrical cloaks, mm 
for elliptic cloaks, [53] 05) for polygonal cloaks, and [541 29] for time domain cloaks). However, as 
commented in mi, such a condition was sufficient but not necessary. Indeed, it was shown in [48J . 
the PEC or PMC could not lead to an independent, meaningful boundary condition for circular 
cylindrical cloaks in polar coordinates. 

Following the spirit of [35], we next introduce the essential “pole” conditions associated with the 
singular transformation (3.2). As we will see, it is important to use the “local” coordinate system 
(r,n) in (3.4) to decompose the differential operators and then carefully study the singularity. 


Proposition 3.1. Let (t,ti) be the “local” coordinate system in (3.4). The CBCs take the form 
X7 t u + = t ■ Vm + = 0 at T+; V„« _ = n ■ X7u~ =0 at (3-12) 

where u + = and u~ = u\qp . 
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Recall that by (2.7), H + = H\ b \qp can be expressed as 


H + = (H+,H+,oy = , - (C 12 u+ + C 22 u +, -C n u+ - C l2 u +, 0)* 


Proof. In order to achieve spectral accuracy in simulations involving singular transformations, we 
follow |3S] (also see [32]) and require that the well-behaved and finite electromagnetic fields in the 
original coordinates r = (£,y) must be well-behaved and finite in the new coordinates r = (x,y). 
We therefore apply this principle to the magnetic field in the cloaking layer, and show that the 
essential “pole” condition of the transformation (3.2 1 takes the form 

V T u+=0 at r*\ (3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 


Then we have 




1t+ = TiVi-U + - T 2 V„M , Uy = T 2 Vt-U + + TiVnlt . 


Inserting (3.6(-(3.8) and (3.15) into (3.14) and collecting the terms, we obtain 

' r — Ri i di?i t 2 \^ + tLRi t± ^ + r 2 ( . ( di?i \ 2 1 > 


U+ (r — K\ dK!T 2 \ + diii Ti + t 2 ( 

=-( K —-r, + —-)X,„u+ + —-V r u + (r+(—) -)V r 


„+ (r-Ri dtflTA 


dl?i T\' 


VtiU+ + ~W~^' rU+ ~ 
dfc 1 r 


Tl 


( /di?r \ 2 1 \ 

i( r+ (^) r) V ^ 


r — R 


By (3.5), 


di?i TiX + r 2 y 

An ’ 

dt/ T2X — T\y 


which is uniformly bounded in f2“ (note: r 2 a; — T\y = 0 implies the side of the polygonal passes 
through the origin, which is not possible). Thus, H + is a finite field in f2“ , if and only if the first 
condition in (3.12) holds. 


Let H = H\^p . From (3.4), (|3.6|)-(|3)8|) and (3.14), we obtain 


n x (H+ H 


I V—Ri 




(0, 0, —n ■ (CVr+) + S7 n u-y\ r=Ri 


— (°’ 
iujUo \ 


n r — R\ , 1 dRi , ' * 

0,-V rl u + H- 

r r ad 


(3.19) 


r—Ri 


where the restriction at r = Ri means that we approach the cloaking boundary from the inside and 
the outside. Suppose that V„m + is finite. Then by ( |3.13| ), 

1 


n x 


(H + - H 


I r—Ri 




(0,0, V T 


I r—Ri ' 


By imposing the second condition in ( 3. 12[ ) , we find 

H~ x n = . (V x E~) x n = 0 at F, 


lLOHo 


(3.20) 


(3.21) 


and by (3.19), the tangential component of the magnetic field H is continuous across the cloaking 
boundary. □ 

Remark 3.2. Weder 


proposed CBCs for point-transformed 3D invisibility cloaks: 
E + x n = H + x n = 0 at dK + ; 

(V x E~) ■ n = (V x H ) • n = 0 at dK_, 


(3.22) 

(3.23) 


where K is the cloaked region. These allowed for the decoupling of the governing equations of the 
inside and outside, and the spherical cloak was considered as a particular application. It is important 
to remark that (3.22)-(3.23) are not applicable to the 2D polygonal cloak, as 

E + x n = (n, 72 , O ) 4 u + 7 ! 0 at Tj. 
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Indeed, E + = (0,0, u + Y does not vanish at the outer cloaking boundary. Moreover, the condition 
(3.23) is different from (3.21). Notably, the CBCs in (3.121 also leads to the decoupling of the inside 
and outside, as we will see below. □ 

With the new CBCs in ( |3.12 ), the governing equations can be decoupled into 

V • (C(r) Vu + (r)) + k 2 n(r)u + (r) =0 in B R \ fF , 

Transmission conditions in (3.11); \/ T u + = 0 at 


and 


d r u + — ^r[m + ] = h at dB R , 


A u + k 2 u =0 in fF ; \/ n u =0 at IF.. 


(3.24) 

(3.25) 

(3.26) 

(3.27) 


Remark 3.3. It is standard to show that (3.27) has a unique solution u (r) = 0, if fc 2 is not an 
eigenvalue of —A in fF with homogeneous Neumann boundary condition. □ 

3.3. Treatment of singularities in spectral-element discretisation. For accurate simulation, 
it is advisable to build the boundary condition V,-it + =0 at r“ in the spectral-element solution 
space. Naturally, we adopt the partition in (3.10), and consider for example f l)p = A p B p BA , where 
A p B p has vertices (aq, yi) and (aa2,2/2)- Suppose that A p B p is mapped to y = —1 via the Gordon-Hall 
transform (2.24|, namely, 


r = xr(£.-l) = ( 


x 2 -xi ^ | aq + x 2 2/2 ~ 2/1 ^ | 2/1 + 2/2 


), £€[- 1 , 1 ], 


which yields 


dxi = 


x 2 - aq 


') ~ 


2/2 - 2/1 


, d x r) = d y tj = 0. 


(3.28) 

(3.29) 


One verifies readily that 
0 = X7 t u + \. 


= {rid x (, + T 2 dy^)d i u e (^, -1 ) + (ti d x r] + r 2 d y y)d v u e (£, -1) 
4 %F(£,-1), 


\/(aq - x 2 ) 2 + (?/i - y 2 ) 2 


where u e = u + (x e (£, —1)). This leads to the corresponding boundary condition in (£, ^-coordinates: 


%F(£,-1) = 0, ^ e [—1,1]. 


Accordingly, we modify the approximation space in (2.19) as 

Vfi = {v £ C(B r \ fF.) :v{r)\ Q e = v(x e (t,y)) € V 2 N , 1 < e < E and 


(3.30) 


(3.31) 


V-ru(r)| n| , nr p = 5 ? u(x e (£, -1)) =0, 1 < e < E T ), 

To meet the boundary condition at FF, we 


3.1 


where E = 12 and Et = 6 for the setting in Figure 
modify the tensorial nodal basis in (2.21): 

^oo=lo(v), Aj = h(O l j(v)> 0 < ? < A, 1 <j<N, (3.32) 

and one verifies readily that 

^Foo = 0, d^ ij \ ri= _ 1 =l'i(t)l j (-l) = 0, 0 <i<N, 1 <j<N. (3.33) 

With such a modification, the singularity can be absorbed by the basis and the spectral-element 
scheme (2.20) can be implemented as usual. 
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3.4. Simulation results for perfect polygonal cloaks. We now provide some numerical re¬ 
sults and compare them with results obtained from the finite-element-based COMSOL Multiphysics 
package. Assume that the incident source is a TE plane wave with an incident angle Oq'. 


i M) = 


Akrcos(O-Oo) _ 


= £ i m J m {kr)e' 


im(6—0o) 


(3.34) 


| m |=0 


so h in (1.9) takes the form 


h = d r u in - &k[u in ] = ifccos(0 - 9 0 )u in - ^ i m J m (kR)T m e lrn8 °, 

| m |=0 


where T m is defined in (1.8). In the following tests, we consider a pentagonal cloak where the polar 
coordinates of ,E are 6 = 7t/5, 37t/7, 87t/9, 67t/5, 77t/4, respectively, and r = 0.7 for all 

vertices. We take p = 0.7 in (3.3) and R = 1.0, and truncate the infinite series in (3.34) with a 
cut-off number M = 60. 




(b) Contour (FEM) 


U 


(c) SEM electric field 




(d) Contour (SEM) 




(e) Profiles of the real and imaginary parts of the electric field along 9 = 0 by SEM 


FIGURE 3.2. A comparison study: SEM versus FEM, where 9 0 = 0 and k = 40. In 
SEM simulation, 40 x 40-grid is used for each element with a total DOF: 16, 000. In FEM 
simulation, a total DOF: 726, 933 is used. 


In Figure 3.2 we plot the electric field distribution (real part) for k = 40 obtained by (i) SEM 
with N = 40 in each element and the degree of freedom (DOF): 16, 000, and (ii) FEM with with a 
total DOF: 726, 933 (in order to obtain reasonable results). Observe that the magnitude of the field 
obtained from FEM is about 1.05 (see the colour bar of (a)), which is expected to be 1 as shown in 
(c). As a result, the wavefront of the field in the contour appears blurred, while that of the SEM 
is very accurate. Indeed, the use of exact DtN boundary condition and new CBCs allows us to 
simulate the ideal cloak very accurately. 

We further challenge SEM with higher wavenumber k = 80 and oblique incident angle 9q = 7 t/4 
(see Figure 3.3). We depict in (a) the electric field distribution (real part) with cut-off number 
M = 80 and N = 80 in each element with the same geometric setting in Figure |ih2| (c). Again, the 
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highly oscillatory oblique incident wave is perfectly steered by the cloaking layer and completely 
shielded from the cloaked region. Apart from plotting the electric-field distributions, we also depict 
the time-averaged Poynting vector (cf. [33]): S = Re{i£ x H*}/ 2, which indicates the directional 
energy flux density. In (b), we depict the associated Poynting vector fields. We find that the waves 
are again steered smoothly around the polygonal cloaked region without reflecting and scattering. 



(a) Electric field 


(b) Poynting vector 


(c) External source 


FIGURE 3.3. SEM with high frequency wave and external source, (a) Real part of the 
electric field distribution and (b) the related Poynting vector, where k = 80, do = n/4. (c) 
Real part of the electric field with external source (3.361. 


We now add an external source, compactly supported in f2 + as the wavemaker, and turn off the 
incident wave. More precisely, we modify (3.24) and (3.26) as 

V • (C(r)Vu + (r)) + k 2 n{r)u + {r) = f{r) in fl + ; d r u + — £^n[u + ] = 0 at Tr. (3.35) 


In practice, we use the Guassian function in Cartesian coordinates: 

tf ^ ( {.x ~ P) 2 + {y ~ k) 2 \ 

f(r) = a exp(^--J, (3.36) 

where a,/3, re, 7 are tuneable constants. To this end, we take a = 100, /3 = —0.41, k = 0.75 and 
7 = 0.04, so the source at Hr is nearly zero. The plot of the electric field distributions in Figure [373] 
(c) is computed from SEM with k = 40, M = 60 and N = 40 in each element. The anti-plane waves 
generated by the source are smoothly bent and the cloak does not produce any scattering. Observe 
that the waves seamlessly pass through the outer artificial boundary without any reflecting. 


3.5. Numerical study of effects of defects, lossy media and dispersive media. We next 
demonstrate that the proposed SEM provides a reliable tool to study the effect of defects and 
sensitivity to the variation of the media within the cloaking layer. Interesting investigation (mostly 
from analytic point of view) has been devoted to the circular and spherical cloaks (see, e.g., mm 
EH0H5U), but the tools appear non-trivial to be extended to the polygonal cloaks. 


3.5.1. Defects in the cloaking layer. We consider the influence of defects to the perfect polygonal 
cloak. As illustrated in Figure [T4| (a)-(b), a rectangular defect with length a and width b is embedded 
into the cloaking layer. We set C = 1 2 and n = 1 within the defect, so the traditional transmission 
condition ( 2 . 10 ) can be imposed at four sides, if the defect is not aligned with the cloaking boundary. 
In Figure 3.4 (a), we depict the electric field distribution with defect a = b = 0.06 obtained by the 
proposed SEM with k = 40, 9q = 0, M = 60 and N = 40. Observe that even with such a small 
defect, the electric field distribution is apparently disturbed, especially for the forward-scattering 
region, and the magnitude increases approximately up to 1.5 (note: it is 1 for the perfect cloak). 
Also notice that in the back-scattering region, the waves appear not significantly affected, so the 
cloaking effects seem still good. In Figure 3.4 (b), we enlarge the defect and set a = 0.24, b = 0.06. 
The field in both the back and forward scattering regions is deteriorated more. 
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FIGURE 3.4. The effects of defects, loss and dispersion on the polygonal invisibility cloak. 
Defects: real part of the electric field distributions with homogeneous defects with sizes (a) 
a = b = 0.06 and (b) a = 0.24, b = 0.06. Loss: real part of the electric field distributions 
with loss tangents (c) tand = 0.01 and (d) tan 5 = 0.05. Dispersion: fix k c = 40, the real 
part of the electric field distributions with wave numbers (e) k = 39 and (f) k = 41. 


3.5.2. Lossy media in the cloaking layer. Similar to the setting in Ell El for the circular and spherical 
cloaks, we replace the media in cloaking layer by a electric-lossy medium (cf. 35]). More precisely, 
the real electric permittivity e in (2.2) is replaced by a complex electric permittivity (1 + itan<5)e, 
where tan S is a tunable constant termed as loss tangent to quantify the absorptive property of the 
medium. Note that this replacement only brings about the modification of the two-dimensional 


Helmholtz equation (2.8) as 

V • 


(C(r) Vu(r)) + k 2 { 1 + itan<5)n(r) u(r) = 0. 


In Figure 3.4 we depict the electric field distributions with loss tangent tan 5 = 0.01, 0.05 in (c) and 
(d), respectively, where we take k = 40, 9q = 0, M = 60 and N = 40 in each element. Observe that 
in the first case, the effect of the loss is almost imperceptible. As we enlarge the loss tangent to 0.05 
in (d), the cloaking effect appears good in the backscattering region but is apparently deteriorated 
in the forward-scattering region, which is inevitable because the lossy medium absorbs the forward- 
travelling wave power. We point out that similar phenomena were observed for the circular and 
spherical cloaks in mm- 


3.5.3. Drude model and dispersive media in the cloaking layer. Based on the form-invariant coordi¬ 
nate transformation, the polygonal cloak can perfectly conceal arbitrary objects inside the interior 
polygonal domain. However, as an ideal cloak, the material parameters are dispersive, and per¬ 
fect invisibility can only be achieved for a single frequency, known as the “cloaking frequency” (cf. 
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[Ml ED E]). It is of much physical relevance to study the response of an ideal cloak to a non- 
monochromatic electromagnetic wave passing through such a dispersive cloak. The investigation 
along this line has been very limited to mostly analytic treatments of circular and spherical cloaks 
(cf. [U[5Tj)- We demonstrate that the proposed SEM offers an accurate means to understand some 
interesting phenomena of a nonmonochromatic wave interacting with a polygonal cloak. 

Following the procedure in [Mjj we start with diagonalizing the symmetric matrices e and fi in 
i.e., 

e = fi = PAP\ A = diag(Ai, A 2 , A 3 ), (3.37) 

where P = (Pij)i<i,j <3 is an orthonormal matrix (with Pj 3 = P 3j = 0 for j = 1, 2, and P 33 = 1), 
and the eigenvalues are 

. O n + O22 + \J {On + C22) 2 — 4 Cn + C 22 — \j (On + C22) 2 — 4 

Ai =-^-> A 2 =--, A 3 = n. (3.38) 


From (3.6)-(3.7), we have 


C n + C22 — 


r — R\ 


1 / 2 fdRi\ 2 \ r — R\ 

+ ^~RCV + \-W) )^— + 


1 — i?i 


> 2 , 


(3.39) 


which implies Ai > 1. However, A 2 and A 3 are less than 1 for some r € (P 1 ,i? 2 ). Based on the 
principle in GHEE], we modify A 2 and A 3 by using the Drude model (cf. [55] 1. More precisely, let 
uj c > 0 be the “cloaking frequency”, and define 

, ,2 


Ai := A i(r,u>) = 1 - 


P,l 


with := w c (w c + iyi)(l - Aj), * = 2,3, 

w(w + i7 i) 


(3.40) 


where {7i}f =2 are given collision frequencies, and {w Pi i}i =2 are known as the plasma frequencies. 
For notational convenience, we define 

w c (w c + by i) 


, so Aj = 1 + Pi (A,; - 1), i = 2,3. 


(3.41) 


w(w + byi) 

Denoting A = diag(Ai, A 2 , A 3 ) with Ai = Ai, we then replace the material parameters e and fi in 


(2.2), respectively, by 


e = ft, = PAP 4 = 


C O 4 

0 fi 


where by a direct calculation, we have 

C n O 12 


C = 


C 12 C, 


22 


and 


= C + (1 — /3 2 )(1 — A 2 ) 
n = 1 + p 3 ( A 3 — 1) = A 3 . 


p2 

-'12 


P12P2 : 

p2 

r 22 


One verifies readily from ( |3.42| ) that 

det(e) = det(ft) = A 3 A 2 A 3 = A 3 A 2 n = det(C) fi. 


Thus, using the fact AiA 2 = 1 (cf. (3.38)), we obtain from (3.41) and (3.45) that 
det(C) = AiA 2 = Ai(/3 2 A 2 + 1 — /3 2 ) = P 2 + (1 — /S 2 ) Ai. 


Accordingly, we find that the counterpart of (2.6) becomes 


ft 1 = e 1 = 


C22 —C 12 1 

-C V2 Cn ' 
0 On 


C, 


where C„ = ?2 + (1 _ a)Ai . 


(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 
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for i,j = 1, 2. Using (2.7) with C. t j and fi in place of C t] and n, we obtain the new model defined in 
the cloaking layer: 

V • (C(r, w) Vu(r)) + k 2 h(r, w) u(r) = 0, (3.48) 

where C = (Cij)i<ij< 2 , and k = ojy/eoPo as before. 


Remark 3.4. Observe from (3.41) that if oj = w c , then /3j = 1 and Aj = Afor i = 1,2. Thus, in this 


case, (3.48) reduces to (|2.8[) and C(r,w c ) is singular at the cloaking boundary r = R±. However, if 


to ^ io c (so P 2 7 ^ 1), then C(r,w) becomes regular at r = i?i. Indeed, by (3.47) 

(r - Ri)Cij 


° ij 02(r-Ri) + (l-02)(r-Ri) Ar' 


In fact, one can verify that if /?2 7 ^ 1, 


lim (r — Ai, A 2 ) all exist. 

!—>Ri 1 ’ 


Thus, we can claim from (3.43) and the above that C(r,u>) is well-defined at r = R\. In view of 
this, the CBCs can not be applied. Here, we follow gUHj and impose a PMC shell instead. □ 

In the computation, we take uj c = k c /y/e 0P0 with k c = 40, and "fi/y/eoPo = 0.0001 for i = 2,3. 
In Figure [3T] (e)-(f), we plot the electric field distributions with k = 39 and k = 41 illuminated by 
plane wave in ( 3.34| ) with incident angle do = 0 and the cut-off number M = 60 and TV = 45 in each 
element. In contrast with Figure 3.2 (c) (where perfect cloaking effect can be obtained for k c = 40), 


we observe from Figure |3.4| that the electric field distributions are affected and distorted in both 
cases (i.e., k = k c zt 1), in particular, more severely when k < k c . Indeed, similar to the phenomena 
observed in EUd for circular and spherical cloaks, the incident wave with frequency slightly deviated 
below fc c , the field after the wave passes the cloak is dissipated and a large shadow appears in the 
forward scattering region. While for the incident wave with frequency slightly deviated above k c , 
the field in most part of the cloaking layer does not change much, except for the part close to the 
cloaking boundary r = R\ 1 and the field behind the cloak is reinforced. This can be regarded as the 
frequency shift effect as in EU- 

4. Accurate simulation of electromagnetic concentrators and rotators 

In this section, we further apply the efficient spectral-element solver to accurately simulate the 
electromagnetic concentrators and rotators. 

4.1. Polygonal concentrators. The electromagnetic concentrator aims at intensifying electro¬ 
magnetic waves in a certain region, which play an important role in the harnessing of light in solar 
cells or similar devices, where high field intensities are needed. 

Here, we are interested in the polygonal concentrator with a configuration similar to the polygonal 


cloak illustrated in Figure 4.1 (b), where EM waves are expected to be concentrated in the interior 
convex polygonal region AU . It is accomplished by a coordinate transformation that maps the 
“polygonal annulus” in Figure 4.1 (a) to the “polygonal annulus” in Figure 4.1 (b), where the 


interior portion of the latter has larger area. More precisely, the polygonal concentrator is mapped 


from the same structure in Figure 4.1 (b) but with a different ratio (see Figure 4.1 (a)): 

OB 0 


P = 


OA 0 
OA OB 


0 < p < p < 1 , 


(4.1) 


where p is defined in (3.3). Then the ratio p/p is known as the rate of concentration. For notational 
convenience, we define 

1 1 ~ P 

Q-= 1 - 1 -=■ 

1 ~ P 

The corresponding coordinate transformation takes the form (see, e.g., [ID]): 


(4.2) 
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(d) Partition of Bj j 


FIGURE 4.1. Schematic geometry of a polygonal concentrator and a circular rotator, 
(a) The polygonal domain in the original coordinates ( x,y ). (b) Through the coordinate 
transformation (4.31, the polygonal domain = A 0 B 0 ---F 0 is compressed into the 
polygonal domain IF that forms the concentration region. Consequently, the original 
polygonal annulus domain \ Cl°_ in (a) is expanded into the polygonal annulus . (c) 
Through the coordinate transformation (4.131, points in the circular annulus a < r < b 
are rotated with a fixed angle 0\. (d) The computational mesh for the circular rotator. 


(i) In 


P 

r = —r, 
P 


f £ [0, Ri], r £ [0, i?i], 
r £ [Ri,R 2 ], r £ 


(4.3) 


J = (l — g)r + qR 2i 
and 9 = 0 £ [0, 2n); 

(ii) In = Br \ the transformation is identity: r = r, 9 = 6. 

Here, R t , i = 1,2, are the same as in ( |3.5[ ) and = p/pRx. Using (2.5), we can derive the 
coefficients C and n as follows (see Appendix |A|): 

(i)a In U0, 

C = I 2 , n = p 2 /p 2 . (4.4) 

(i )b In Hd, 


C n = 


r — qR 2 x 2 


g f g dR 2 x 2 2 xy \ dR 2 


+ 


r — qR 2 y 2 
U'22 —-y + 


r — qR 2 \ r d9 r 2 r 2 J d9 r — qR 2 r 2 
g f gdR 2 y 2 2xy\dR 2 r x 2 


r r“ r — gR 2 \ r d9 r 2 r 2 J d9 r — gR 2 r 2 ’ 


C12 — 


r - gR 2 


xy 

o ' 


Q 


g dR 2 xy y 2 \ d R 2 


and 


r - gR 2 Jr 2 r - gR 2 \ r 


r - gR 2 


r d 9 r 2 


d 9 


r(l — g) 2 

(ii) In we have C = 1 2 and n = 1. 

In summary, the governing equation for the polygonal concentrator reads 
V • (C(r) Vw(r)) + fc 2 n(r)tt(r) = 0 in Br, 

[u] = |CVu]=0 at PUT 0 , 

d r u — 5 r[u] = h at 8Br. 


(4.5) 

(4.6) 

(4.7) 

(4.8) 


(4.9) 

(4.10) 

(4.11) 
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Note that in the interior polygon, ( |4.9[ ) becomes the Helmholtz equation: 

in il T> _, 


o 2 

A u + ~^k 2 u = 0 
P 


(4.12) 


where the ratio p 2 / p 2 > 1 is a constant. Therefore, the coordinate transformation enlarges the 
wavenumber k that produces the effect of concentration. 

We can implement the spectral-element solver based on the partition of the computational domain 
as with the previous application. However, different from the previous case, the interior region is part 
of the computational domain, where a normal transmission condition is imposed along its boundary 
(see (4.10)). Below, we provide some numerical results with the setting: a square concentrator 
centred at the origin with length of each side 1.2 and the parameters: p = 1/3, p = 2/3 in (3.3) and 
(4.1) and R = 1.0. We set the cut-off number M — 60 in the DtN operator. In Figure [O] we depict 


the electric field distributions and the associated time averaged Poynting vectors illuminated with 
different incident angles ((a)-(b): do = 0 and (c)-(d): 9q = 7 t/4 ) with k = 40 and the grid N = 50 in 
each element. It can be seen that the electric field and energy flux are smoothly concentrated into 
the inner concentration region , and the field outside is not affected regardless of the incident 
angle on the concentrator. 



FIGURE 4.2. The real part of the electric field distributions and associated Poynting 
vectors for square concentrators with (a)-(b): do = 0 and (c)-(d): do = 7r/4, respectively. 


4.2. Circular rotators. In contrast with the invisibility cloak and concentrator, the electromag¬ 
netic rotator is based upon a coordinate transformation of the angular variable rather than the 
radial variable in (3.1). 


As illustrated in Figure 4.1 (c), the domain is a disk of radius r = b, which encloses a 
concentric disk of radius r = a < b. The waves are expected to rotate with a fixed angle in the 
interior disk. This can be realised by the coordinate transformation (cf. [5]): 


r = r, 6 = 9 + 9\, 

„ n x s (b) - s(r) 
r = f, 9 = 9 -{—/A- T\9\i 


0 < f < a, 9, 9 € [0, 2tt), 
a < r < b, 9, 9 £ [0,27r), 


(4.13) 


s(b) — s(a) 

where s is any smooth function such that s(b) ^ s(a). As before, the transformation is identity 
exterior to . 

Define 

s(r)' 


K = 


■0i. 


s(b) — s(a) 

Working out the material parameters as before (see Appendix |A|, we have 
1 r 2 + 2 nxy + n 2 y 2 —nx 2 — n 2 xy + ny 2 
—kx 2 — n 2 xy + ny 2 r 2 — 2nxy + k 2 x 2 




n = 1 for a < r < b, 


(4.14) 


(4.15) 
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and 


C = I 2 , n= 1 for 0 < r < a, b < r < R. 


(4.16) 


In Figure 4.1 (d), we illustrate a partition of the computational domain Br. Together with the 


standard transmission conditions in (2.12), we can implement the spectral element scheme as the 


previous cases with a similar partition of the computational domain (see Figure 4.1 (d)). In the 
computation, we set a = 0.3, b = 0.7 and s(r) = r in (4.13) and choose R = 1.0, M = 60 and 


N = 40 in each element. In Figure |4~3} we fix k = 40, and plot the electric field distribution (real 
part) and the corresponding time averaged Poynting vector with the same incident angle 9q = 0 
and different rotation angles ((a)-(b): 9\ = 7r/4, (c)-(d): 9\ = 37t/ 4). We find that the electric field 
distribution rotates its direction by 7r/4 in Figure [473] (a) and the power flux (b) flows with the same 


direction in the closed region r < a. It can be observed in Figure 4.3 (c)-(d) that even for a very 
sharp rotation angle 9\ = 37r/4, the field rotates exactly by 37 t/ 4 angle without introducing any 
scattering wave outside. 



FIGURE 4.3. The real part of the electric field distributions and associated Poynting 
vectors for circular rotators with (a)-(b): 9i = 7t/4 and (c)-(d): 9\ = 37t/4, respectively. 


Concluding remarks 


In this paper, we presented an accurate and efficient spectral-element solver for time-harmonic 
Helmholtz equations in general inhomogeneous and anisotropic media. We focused on several appli¬ 
cations arisen from transformation electromagnetics which included the polygonal invisibility cloaks 
and concentrators, and circular rotators. We introduced new ideas of how to seamlessly integrate 
local elements and global DtN boundary condition. This can shed light on the three-dimensional 
simulation where the DtN BC involves spherical harmonic expansions. We proposed new cloaking 
boundary conditions for accurate simulation of perfect polygonal invisibility cloaks. The proposed 
method also provided a reliable alternative to the analytic tools to study the interesting phenomena 
when defects and other media were embedded or placed in a perfect cloak. 


Appendix A. Derivation of material parameters for polygonal cylindrical cloaks, 

CONCENTRATORS AND CIRCULAR ROTATORS 


As illustrated by (3.1) in Section [ 3 J the transformation is given from polar coordinates (r,9) 
(of the original Cartesian coordinates (x,y)) to polar coordinates (r, 9) (of the physical Cartesian 
coordinates ( x,y )), so by the chain rule, the Jacobian matrix J cn can be computed by 


d$x 

djjX 


d r x 

d 9 x 

d?r 

d e r 

OxX d$r 

dxV 

dyy 


d r y 

dgy_ 

d f 9 

-1 

px9 dtjO 
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It is clear that 


and 


Denote 


d r x 

dgx' 

cos 9 

—r sin# 

det ^ 

d r x 


pry 

dgy_ 

sin 9 

r cos 9 ’ 

p r y 

II 

■ 

dxr 

dp 

cos 9 

sin0 

det | 

' 'dxr 

<V] \ _ 1 

dp 

dy9_ 

— sin 9/r 

cos 9/r\ 

v dp 

dy9\ J f 


d?r dgr 

d?e 8 § e\ ’ 


(A.2) 

(A.3) 

(A.4) 


which is determined by specific coordinate transformations. With a direct calculation, we derive 
from (2.5) the material parameters C and n for the general transformation in (3.1): 


Cu = r 3 d( ! t( -j) (( xd ? r - r ydf8) 2 + ^ (xd?6 - ryd g 6) 2 ^j , 
C22 = r 3 de ? t (j) {(y dfr + rxd ?0) 2 + ^ {ydfr + rxdgO) 2 ^ , 


(A.5) 

(A.6) 


Cy> = 


r 3 det(J) 


(*v($r + - r 2 d?d ^dp) + ^ / )r (v*d f rd f 6 + dgrdgO)), (A.7) 


and 


det(J cn ) rdet(J) 


For the polygonal cylindrical cloak, we derive from (3.2) that in f2“ , 

1 


r = 


1 ~ P 


(r — R±), e = e, 


which leads to 


J = 


1 - p d e R\ 
0 1 


, det( J) = 1 — p. 


(A.8) 

(A.9) 

(A.10) 


Recall that dgRi can be worked out by (3.5). Thus, the material parameters C and n in (3.6)-(3.9) 
can be obtained by substituting (I A.9h- (|A. 10h into ( A.5h-(|A.8 ). 


Similarly, for the polygonal cylindrical concentrator, we obtain from the transformation (4.3) that 
in fii, 

-(r - qR 2 ), 9 = 9 , 


and 


J = 


1 ~ Q 

1 - B gd e R 2 
0 1 


, det(J) = l — g. 


Inserting (A.11) and (A.12) into (A.5)-(A.8), we obtain C and n in (4.5)-(4.8). 
For the circular rotators, we have from (4.13) that 

a a s ( b )~ s ( r ) Q 
r = r, 9 = 9 -—- —9 1 , 


s(b) — s(a ) 


and 


~d?r 

d 6 r 


' 1 O' 

d?9 

1 


— K 1 


(A.ll) 
(A.12) 

(A.13) 
(A. 14) 


where k is defined in (4.141. Then we can compute the material parameters (4.14) in a similar 
fashion. 
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